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ABSTRACT 


I 

|  The  purpose  of  this  investigation  is  to  determine  the  in-plane  thermal  stresses  for  a 

symmetrically  laminated,  50”x  1 2”x0. 1 9”composite  plate  with  temperature  dependent 
material  properties.  For  this  study  only  in-plane  stresses  are  investigated.  The  in-plane 
equations  of  motion  are  solved  exactly  using  a  stress  function,  and  the  resulting  compati¬ 
bility  equation  is  solved  approximately  using  the  Galerkin  method.  This  investigation  also 
serves  as  proof  of  concept  for  the  variational  method .  This  method  produces  accurate 
results  while  being  less  rigorous  in  a  computational  sense  than  the  high  degree  of  free¬ 
dom  finite  element  model  required  to  solve  the  same  problem.  Variations  with  lamina  ori¬ 
entation  and  multiple  layered  laminates  are  investigated.  Results  are  given  in  terms  of  the 
in-plane  force  resultant.  The  baseline  case  f  or  this  study  was  aluminum  and  the  in-plane 
force  resultant  at  the  center  of  the  plate  was  calculated  to  be  60.251  lb/in.  The  exact  solu¬ 
tion  for  the  in-plane  force  resultant  at  the  center  of  the  plate  is  59.667  Ib/in,  a  difference  of 
less  than  1  percent.  Based  on  these  results  additional  investigations  were  accomplished 
for  composite  plates.  The  results  from  this  study  will  be  used  as  parametric  data  by 
NASA-Dryden  in  verifying  finite-element  codes  and  will  aid  in  experimental  analysis  of 
thermal  loading  on  composite  plates.  Furthermore,  this  study  should  serve  as  a  reference 
to  an  investigation  that  considers  thermally  stressed  laminates  with  bending  and  exten- 
sional  coupling. 
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1.  INTRODUCTION 

1.1  MOTIVATION 

Hypersonic  atmospheric  vehicle  research  is  again  becoming  an  area  of  great  interest. 
The  push  to  develop  a  vehicle  which  can  take-off,  land,  and  be  operated  in  a  fashion  simi¬ 
lar  to  current  aircraft  is  one  of  the  interests  in  this  area.  Current  concepts  incorporate  com¬ 
posite  materials  (possibly  ceramic  and/or  metal  matrix)  as  one  of  the  principal  structural 
components.  Studies  at  NASA-Dryden  are  being  conducted  to  verify  the  capabilities  of 
composite  materials  in  such  an  application.  The  results  of  this  study  and  others  in  this  area 
will  be  used  by  NASA-Dryden  as  an  analytical  tool  to  verify  finite  element  methods  for 
heated  structures  and  to  aid  in  experimental  investigations. 

In  addition  to  the  mechanical  stresses  present  in  composite  materials,  the  thermal 
stresses  must  also  be  investigated,  especially  for  hypersonic  flight  applications.  For  com¬ 
posite  plates,  coupling  of  the  extensional  and  bending  deformations  is  usually  present 
However,  for  this  study  only  extensional  deformations  are  investigated.  The  purpose  of 
this  study  is  to  develop  a  mathematical  model  for  the  in-plane  thermal  stresses  of  a  sym¬ 
metrically  laminated,  rectangular,  composite  plate  with  stress  free  boundary  conditions. 

For  the  present  study,  the  classical  Kirchoff  thin-plate  theory  is  used  and  the  in-plane 
equation  of  motion  is  expressed  in  terms  of  a  stress  function.  The  resulting  partial  differ¬ 
ential  equation  is  solved  using  Galerkin’s  method.  The  temperature  distribution  for  this 
study  is  symmetrical  and  the  relation  of  the  material  properties  to  temperature  is  assumed 
to  be  linear. 

1 .2  THEORETICAL  BACKGROUND 

1.2,1  Stress-Strain  Relations  for  Composite  Plates 

Laminated,  fibrous  composite  plates  are  typically  made  of  stacked  layers  of  fibers, 
each  layer  having  all  its  fibers  aligned  in  the  same  direction.  However,  the  alignment  most 
often  varies  with  each  layer.  In  order  to  characterize  the  stiffness  of  the  resulting  plate,  it 
is  important  to  be  able  to  write  the  stress-strain  relation  for  each  layer.  In  particular  the 
stress-strain  relationship  must  be  known  with  respect  to  an  axis  system  which  is  not  neces¬ 
sarily  aligned  in  the  same  direction  as  any  of  the  fibers  in  any  given  layer. 

The  generalized  Hooke’s  Law  relating  stresses  to  strains  can  be  written  as 


-  cijej 


i,j  =  1...6 


(EQl) 


where  Cjj  is  the  stiffness  matrix  and  the  O-,  are  the  stress  components  and,  e j  are  the  strain 
components.  The  stiffness  matrix  has  36  constants,  however,  due  to  symmetry  the  stiff¬ 
ness  matrix  is  populated  by  21  independent  constants  (REF  2).  Thus,  in  simplest  terms  the 
stress-strain  relations  for  any  linear  elastic  material  is  given  by  the  general  expression 
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The  relations  in  (EQ  2)  are  referred  to  as  characterizing  anisotropic  materials  since 
there  are  no  planes  of  symmetry  for  the  material  properties.  An  example  of  this  is  a  thick 
laminated  plate  (REF  1). 

If  there  are  three  orthogonal  (mutually  perpendicular)  planes  of  symmetry,  such  as  for 
a  single  layer  of  Fibers,  the  material  is  said  to  be  orthotropic.  The  stress-strain  relation  for 
an  orthotropic  material  is  similar  to  (EQ  2),  however,  now  there  are  only  9  independent 
constants  instead  of  21.  The  remaining  constants  are  zero.  The  nonzero  constants  are  Cj  i» 
C12,  C13.  C22*  C23,  C33,  C44,  C55,  C55.  Furthermore,  in  work  done  by  Lempriere  it  was 
shown  that  the  stiffness  matrix  must  be  positive-definite  (Cjj>0  for  Cl}  0).  This  condi¬ 
tion  is  a  consequence  of  the  requirement  that  the  sum  of  the  work  done  by  all  stress  com¬ 
ponents  must  be  positive  in  order  to  avoid  the  creation  of  energy  (REF  2). 

1.2.2  Two  Dimensional  State  of  Stress 

Two  dimensional  stress  or  plane  stress  is  an  area  of  great  interest  in  many  different 
fields.  One  concept  of  plane  stress  is  a  flat  plate.  Flat  plates  are  used  as  a  model  in  a  vari¬ 
ety  of  analytical  approximations  for  otherwise  untractable  problems.  In  many  problems, 
choosing  the  z-direction  to  be  normal  to  the  plate  surface,  the  stresses  oz,  xyz,  and  x^  are 
small  with  respect  to  the  other  stresses.  Thus,  the  stress-strain  relationship  for  an  orthotro¬ 
pic  material  in  a  state  of  plane  stress  is: 


°x 
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(EQ3) 


This  is  the  relationship  relevant  to  a  layer  of  fibers  in  a  composite  laminate,  the  fibers 
being  aligned  in  the  x-direction  (REF  1). 

Instead  of  considering  the  fibers  to  be  oriented  along  the  (x,y)  axis  system,  let  us  con¬ 
sider  the  principal  axes  of  a  given  lamina,  with  the  1 -direction  parallel  to  the  fibers  and  the 
2-direction  perpendicular  to  the  fibers.  In  doing  this  the  C’s  in  (EQ  3)  are  replaced  with 
Q’s  and  are  now  referred  to  as  reduced  stiffnesses.  Furthermore,  x  and  y  are  replaced  with 
one  and  two  respectively.  In  terms  of  engineering  constants  the  reduced  stiffnesses  are: 
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(EQ4) 
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Where  Ej  is  the  modulus  in  the  1-direction,  Vj2  is  the  Poisson’s  ratio  for  contraction  in  the 
2-direction  due  to  extension  in  the  1 -direction,  and  G12  is  the  shear  modulus  relating  shear 
stress  in  the  2-direction  to  shear  strain  in  the  1 -direction. 


The  preceding  stress-strain  relations  are  the  basis  for  the  stiffness  and  stress  analysis 
of  an  individual  lamina  subjected  to  forces  in  its  own  plane.  These  relations  are  essential 
in  the  analysis  of  laminates  (REF  2). 

1.2.3  Stress-Strain  Relations  for  a  Lamina  of  Arbitrary  Orientation 


In  most  problems  it  is  necessary  to  solve  for  the  stress-strain  relationship  in  a  direction 
other  than  the  principal  material  coordinates.  An  example  is  a  laminated  plate  consisting 
of  lamina  oriented  at  different  angles  relative  to  the  plate  axis  system.  The  stress-strain 
relationship  for  a  lamina  in  the  principal,  (1,2),  system  transformed  to  a  global  (x,y)  sys¬ 
tem  is  given  by: 


cos20  sir20  ~2sin0cos0 

ai~ 

= 

2  2 

sin  0  cos  0  2sin0cos0 

°2 

2  2 

sin0cos0 -sin0cos0  cos  0-sin  0 

J\2 

(HQ  5) 


where  0  is  the  angle  from  the  x-axis  to  the  1-axis  (REF  2).  In  some  treatments  of  the  sub¬ 
ject  the  quantity  eXy=l/2yXy  is  introduced  and  is  termed  the  tensor  expression  for  shear 
strain  (REF  1).  In  compact  form  the  stress-strain  relationship  for  a  lamina  oriented  at  any 
angle  is: 


(EQ6) 


where  Q^: is  the  transformed  reduced  stiffnesses  instead  of  the  reduced  stiffnesses  Qjj,  and 
k  is  the  k®  layer  of  a  multilayered  laminate.  The  equations  for  QtJ  are  given  in  References 
1,  and  2.  Also  note  that  the  transformed  reduced  stiffness  matrix  has  terms  in  all  nine 
positions  instead  of  the  zeros  present  in  the  stiffness  matrix 


1.2.4  Mechanical  Properties  of  Laminated  Composite  Plates 

In  the  preceding  section  the  stress-strain  relations  for  a  lamina  of  an  orthotropic  mate¬ 
rial  under  plane  stress  oriented  at  any  angle  were  developed.  These  relations  are  useful 
when  dealing  with  laminated  plates  because  of  the  arbitrary  orientation  of  the  lamina. 
Equation  6  can  be  thought  of  as  the  stress-strain  relations  of  the  k®  layer  of  a  multilayered 
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laminate.  The  next  step  is  to  develop  the  stress  and  strain  variations  through  the  thickness 
of  a  laminate. 


1.2.5  Strains  for  a  Laminated  Plate 

The  laminate  is  presumed  to  consist  of  perfectly  bonded  lamina  which  are  not  allowed 
to  slip  relative  to  another.  Furthermore,  the  bonds  are  presumed  to  be  very  thin  as  well  as 
non-shear  deformable.  Additional  assumptions  of  the  behavior  of  the  laminate  are  given 
by  the  Kirchhoff  theory  of  plates  (REF  2). 


The  development  of  the  strains  for  a  laminated  plate  are  given  in  Reference  2.  Ulti¬ 
mately,  the  laminate  strains  are  reduced  to  ex,  £y,  and  yxy,  (i.e.Yxz-^Yyz-^^O)  by  virtue  of 
the  Kirchhoff-Love  hypothesis.  Using  the  derived  displacements  u  and  v  for  the  x  and  y 
directions  respectively,  the  strains  in  matrix  form  are 


e  = 
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(EQ7) 


Here  e°x,  e®  ^xy,  are  the  “midplane  strains”,  Kx  and  Ky  are  the  “curvatures”  and  Kxy  is 
the  “twist”  (REF  1).  Now  applying  (EQ  6)  and  (EQ  7)  one  can  determine  the  stresses  in 
the  k111  layer  in  terms  of  the  midplane  strains,  curvatures,  and  twist  using  the  relation: 
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In  this  problem  only  in-plane  stresses  are  investigated  so  the  curvatures  and  twists  are 
neglected. 


1.2.6  Resultant  Laminate  Forces  and  Moments 

Typically  a  designer  deals  with  forces  and  moments  instead  of  stresses.  The  resultant 
in-plane  forces  are  calculated  by  integrating  the  stresses  over  the  thickness  of  the  plate. 
The  common  load  parameters  for  a  plate  are  the  edge  force  intensities,  Nx,  Ny,  and  Nxy  in 
units  (lb/in),  which  are  obtained  using  the  following  equation: 
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Here  N  is  the  total  number  of  layers,  h  is  the  total  plate  thickness  and  k  is  the  k1*1  layer. 
Since  the  transformed  reduced  stiffness  matrix  in  (EQ  8)  is  constant  for  each  lamina,  upon 
substitution  into  (EQ  9),  it  can  be  pulled  outside  the  integral,  however,  it  must  remain 
within  the  summation  of  force  resultants  for  each  layer.  In  addition,  the  middle  surface 
strains  are  not  dependent  on  z  so  they  can  be  removed  from  under  the  integral  and  summa¬ 
tion  signs  (REF  2).  If  one  defines  a  matrix 
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H  =  I 

*=i 


(EQ  10) 


then  the  in-plane  loads  for  a  laminated  plate  can  be  expressed  as 


M  =  M 


(HQ  1!) 


where  [A]  is  referred  to  as  the  “extensional  stiffness”  matrix.  Since  only  in-plane  force 
resultants  are  evaluated:  |g0  =  (REF  2). 


1.2.7  Stress-Strain  Temperature  Relations 

As  stated  in  the  introduction,  the  purpose  of  this  project  was  to  determine  the  thermal 
stresses  in  a  laminated  composite  plate.  The  equations  developed  to  this  point  include 
only  the  mechanical  stresses  and  do  not  account  for  the  thermal  stresses.  Thermal  stresses 
arise  from  the  expansion  and  contraction  of  fibers  in  a  solid  material.  “For  the  solid  to 
remain  continuous,  a  system  of  thermal  strains  and  corresponding  thermal  stresses  may  be 
induced,  depending  on  the  characteristics  of  the  solid  and  its  temperature  distribution” 
(REF  4). 

For  plane  stress  of  an  orthotropic  material  in  principal  material  coordinates  the  stress- 
strain-temperature  relation  in  terms  of  the  stiffness  matrix  is  given  by:  (REF  2) 
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Here  a  is  the  coefficient  of  thermal  expansion  in  the  principal  material  directions  and  AT 
is  the  relative  change  ir.  temperature  from  n  datum  temperature.  For  an  arbitrary  orienta¬ 
tion  the  coordinate  transformation  accomplished  in  (EQ  5)  must  be  performed.  After  the 
transformation  to  laminate  coordinates  a  thermal  shear  term,  axy,  is  present  where  only  the 
extensional  strains  were  present  before: 
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Substituting  (EQ  13)  into  (EQ  9)  and  defining  the  following  thermal  load  vector 
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one  can  derive  the  following  in-plane  resultant  forces: 
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(EQ  15) 


The  thermal  forces,  NT,  are  true  thermal  forces  only  when  the  total  strains  and  curvatures 
are  zero  (REF  2). 

Rearranging  (EQ  15)  to  solve  for  the  strains  results  in  the  following  contracted  matrix 
equation 


(EQ  16) 


Here  [a]  is  the  inverse  of  the  extensional  stiffness  matrix  [A]. 

This  is  a  relatively  short  description  of  the  classical  laminated  plate  theory  required  for 
the  present  study. 
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The  purpose  of  this  section  to  take  a  step  by  step  approach  to  the  solution,  for  in-plane 
loads,  of  a  thermally  stressed  composite  plate. 

2.1  FORCE  EQUILIBRIUM  EQUATIONS 

“To  obtain  the  exact  stress  field  in  any  given  solid  under  the  action  of  external  loads, 
the  equations  of  equilibrium,  compatibility  equations,  and  the  boundary  conditions  must 
all  be  satisfied”  (REF  4).  Therefore,  in  developing  the  equation  for  thermal  stresses  it  is 
wise  to  siart  with  the  equilibrium  equations. 

The  force  equilibrium  equations  for  a  2-dimensional  state  of  plane  stress  with  zero 
body  forces  are: 


dhJ3> 

dx  dy 


0 


Wx, 

dx 


0 


where  Nx,  Ny,  and  Nxy,  are  in-plane  loads  with  units  of  force/length. 
The  “Airy”  stress  function  F  defined  such  that: 


(HQ  17) 


&  djF  J2F 
dyrdxr  dxdy_ 


identically  satisfies  the  equilibrium  equations. 


2.2  COMPATIBILITY  EQUATION 


(HQ  18) 


The  compatibility  equation  is  an  equation  that  gives  a  relationship  between  the  deriva¬ 
tives  of  the  strain  components  (REF  5).  The  compatibility  equation  for  a  two  dimensional 
strain  field  can  be  derived  from  the  strain  displacement  relations  given  in  section  1.2.5. 
Differentiate  the  third  equation  in  (EQ  7)  with  respect  to  x  and  y  to  obtain 


d\y  =  d3u  (  83v 
dxdy  dy2dx  dx2y 


(EQ  19) 


Substituting  the  first  two  equations  from  (EQ  7)  into  (EQ  19)  one  can  obtain  the  compati¬ 
bility  equation  for  a  two  dimensional  strain  field: 


dxdy  dy2  dx2 


(EQ  20) 


Substitution  of  (EQ  16  &18)  into  (EQ  20)  results  in  the  partial  differential  equation  form 
of  the  compatibility  equation: 


By 


;(* 


Udy2+a'2dx 2 


d2F 


^ ,  a2f  b2f  ,  b2f  a^p  ^ 

j B?lai2Bf  22B?  °26^) 

(ai6a/+a26a?_fl66^y)  =  -^(anN*+anN*+ai<>N*y^ 


(HQ  21) 


BxBy{  16ay 

-^(*l2»l  +  *22HTy  +^26<)  +  ^(*16*1  +  +  *66^) 

2.3  VARIATIONAL  METHOD 


The  phrase  ‘variational  methods’  refers  to  methods  that  make  use  of  variational  princi¬ 
ples  to  determine  approximate  solutions  to  partial  differential  equations.  In  this  investiga¬ 
tion,  a  variational  formulation,  which  will  be  recognized  as  the  principle  of  minimum 
complementary  potential  energy,  is  used  to  approximate  a  solution  to  the  compatibility 
equation  described  in  section  2.2. 


2.3.1  Variational  Formulation 

The  first  step  in  the  variatonal  formulation  is  to  multiply  (EQ  20)  with  all  the  terms  on 
one  side  of  the  equality  with  the  variation  5F  and  integrate  over  the  domain: 

I5F<e*.,,+e,.«-W‘M  =  °  «Q22> 

A 

where  it  is  understood  that  all  the  strains  are  expressed  in  terms  of  the  stress  function  F. 

Equation  22  can  be  written  in  an  alternative  form  by  applying  the  Green-Gauss  theo¬ 
rem  twice  to  trade  the  differentiation  between  F  and  8F  so  they  both  have  the  same  order 
derivatives  (REF  5).  In  applying  the  Green-Gauss  theorem  the  first  step  is  to  rewrite  (EQ 
22)  in  the  form: 

»AdA  = jn*AdS  (EQ23) 

A  S 

where 


k  =  SF(£,.-kI)i+5F(EM-!v,)5 


n  =  nj  + 


ny 


Applying  the  Green-Gauss  theorem  twice  results  in 


(EQ  24) 


f(SF  e  +8F  e  -8F  y  )dA  =  0 

J  K  \xx  y  ,yy  X  *ylxyJ 
A 

Expanding  (EQ  25)  for  the  strains  results  in  the  equation: 

J  [8 F  ( anF  yy  +  -  a16F  ^  +  anNT x  +  al2NTy  +  al6NTxy) 

A 

-4-  SF^  {ax2Fyy  +  an -  o^F^  +  al2N^+  a22Ny  + 

+  SFsy(a\6F,yy  +  a26F^-a66Fjy  +  al6NTx  +  a2ftNTy  +  a66Nly)  ]dA 

2.3.2  Galerkin’s  Method 

The  integral  equation,  (EQ  25),  is  the  variational  equivalent  of  the  in-plane  compati¬ 
bility  equation,  (EQ  20),  and  essentially  represents  a  Galerkin  solution  for  the  partial  dif¬ 
ferential  equation  shown  explicitly  in  (EQ  21). 

The  first  step  in  the  Galerkin  method  is  to  approximate  the  unknown  function  with  a 
suitable  set  of  admissible  functions.  For  the  present  formulation  the  stress  function  F  is 
approximated  as: 

n 

F«  ]£  Ffc/CV  (EQ  27) 

J u=  1 

where  £k  and  T]1  are  required  to  satisfy  the  following  three  conditions: 

a)  be  continuous,  as  required  by  the  variational  principle  being  used 

b)  satisfy  the  specified  essential  boundary  conditions 

c)  be  linearly  independent  and  complete  (REF  5) 

For  this  problem  the  essential  boundary  conditions  are  zero  and  a  series  of  polynomials 
were  chosen  as  the  approximating  function  because  polynomials  can  satisfy  arbitrary 
boundary  conditions,  and  Fy  are  the  coordinates  or  components  of  the  approximation. 

The  assumed  solution  in  (EQ  27)  is  in  the  form  of  a  finite  linear  combination  of  unde¬ 
termined  parameters.  Approximating  the  continuous  function  in  (EQ  25)  by  a  finite  linear 
combination  of  functions  introduces  some  error.  Therefore,  the  solution  obtained  is  an 
approximation  of  the  true  solution  for  the  equation  of  motion  described  in  (EQ  25  &  26). 

The  variation  of  (EQ  27)  is  given  by 

n 

8F  =  £  SF;>CV  (EQ  28) 

*.7=  1 


(EQ  25) 


(EQ26) 
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2.3.3  Matrix  Form  of  the  Compatibility  Equation 

Taking  the  variation  of  (EQ  18)  it  is  clear  that  {8Nx,8Ny,8Nxy}=(SF,yy,8F,xx,-8F,xy}. 
Thus,  the  matrix  form  of  (EQ  25)  is 

J{8N}t{£}<M  =  0  (EQ  29) 

A 

where  {8NT}  is  the  transpose  of  the  matrix  in  (EQ  18).  Substituting  (EQ  16)  into  (EQ  29) 
to  obtain  the  matrix  equation 


|[{8N}r[fl]  {N}  +  {8N}  T[a]  {NT})dA  =  0  (EQ30) 

A 

Since  (EQ  30)  is  in  matrix  form  it  is  necessary  to  put  (EQ  28)  into  matrix  form  also. 
This  is  accomplished  by  using  a  row  vector  of  length  N,  where  N=(n+l)*(n+l)  and  n  is 
the  order  of  the  polynomial  in  (EQ  27)  and  (EQ  28),  and  a  column  vector  of  length  N  in 
the  following  fashion 


f=  XfyCV-  n.n.n2.  Cn2.  ...?V.C2.  -CV1 

*1^00’  pou  ^02«  •••>  Fo*  Fl0,  Fjj,  Fj2 . FIn,  F20,  --F ^ ] T 


(EQ  31) 


Here  the  row  vector  can  be  denoted  by  [H]  and  the  column  vector  by  (F).  However,  since 
(8N)  is  composed  of  second  order  derivatives,  and  it  relates  directly  to  (8F) ,  it  is  neces¬ 
sary  to  differentiate  F  in  (EQ  31)  according  to  the  relation  in  (EQ  25).  Only  [H]  needs  to 
be  differentiated  since  {F}  does  not  depend  on  x  and  y,  therefore 


f„=  m„{F) 

(EQ  32) 

w^iF} 

F,y= 

If  a  matrix,  [B],  is  defined  to  contain  the  three  row  vectors,  [H]  yy,  [H]  xx,  and  [H]  xy, 
then  [B]  would  be  a  3xN  matrix,  and  {N)=[B]{F).  Similarly  {SN}=[B]{8F}.  From  this 
relation  it  is  clear  that  {8N)T={8F)T[B]  .  If  these  results  for  {N}and  {8N}  are  substi¬ 
tuted  in  (EQ  29)  the  matrix  equation  is  given  by 

f({8F}T[B]T[aJ  [B]  {F}  +  {SF}  T [B]  T[a]  { NT})dA  =  0  (EQ33) 

A 

Since  {5F}^is  common  in  both  factors  of  the  integrand  it  can  be  factored  to  one  side.  In 
addition,  for  arbitrary  and  linearly  independent  SFy  (i,j=l,n)  it  follows  that  the  remaining 
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integrand  must  be  equal  to  zero.  If  the  thermal  terms  are  brought  to  the  right  hand  side  of 
the  equation  and  the  stress  function  component  matrix  is  brought  outside  the  integral 
because  it  is  constant  then 

J([£jr[0]  [B})dA{F}  =  —  J  (  [5]  T  [d\  {Nt}  )dA  {3Q34) 

A  A 

Thus  we  obtain  N  linearly  independent  simultaneous  equations  for  the  N  unknowns  in 
{F}.  Once  the  coefficients  {F}  are  determined  the  in-plane  force  resultants  can  be  deter¬ 
mined  at  any  point  within  the  plate  by  using  the  relation  {N}=[B]{F}. 

2.4  NUMERICAL  INTEGRATION 

In  order  to  evaluate  the  [A]  { x } = { b }  problem  shown  in  (EQ  34)  it  is  necessary  to  inte¬ 
grate  the  [A]  and  (b)  matrices.  The  integration  could  be  done  explicitly,  however,  an 
alternate  method  that  yields  exact  results  for  polynomial  functions  is  Gaussian  quadrature. 

Most  numerical  integration  schemes  are  based  on  predetermined  x-values.  However, 
Gauss  observed  that  if  we  remove  the  requirement  that  the  function  be  evaluated  at  prede¬ 
termined  x-values,  then  a  three  term  formula  will  contain  six  parameters:  the  three 
unknown  x-values  and  the  three  weights  (REF  7).  This  should  correspond  to  an  interpo¬ 
lating  polynomial  of  degree  5.  In  addition,  the  tabulated  values  for  the  Gaussian  quadra¬ 
ture  procedure  are  derived  for  a  symmetric  interval  of  integration  from  -1  to  1.  Clearly, 
Gaussian  quadrature  formulas  can  only  be  applied  when  the  function  is  explicitly  known. 
Since  the  integrand  in  (EQ  34)  is  approximated  by  a  polynomial  it  can  be  calculated  at  any 
point  and  is  readily  suited  to  this  method  of  numerical  integration. 

For  this  problem  the  limits  of  integration  were  not  from  -1  to  1.  The  [B]  matrices  are 
in  terms  of  £  and  rj,  with  the  origin  of  this  axis  system  placed  at  the  center  of  the  plate. 
The  limits  of  integration  for  this  axis  system  were  from  -a/2  to  a/2,  and  -b/2  to  b/2,  where 
a  and  b  are  the  length  and  width  of  the  plate  respectively.  Therefore  £  and  r|  were  set 
equal  to  2x/a  and  2y/a.  Now  the  limits  of  integration  are  from  -1  to  1  and  a  factor placed 
in  front  of  the  integral  accounts  for  the  change  of  variable  effects  for  dA=dxdy=^~  dt,dr[ . 


it 


The  purpose  of  this  section  is  to  explain  how  the  program  to  calculate  the  coefficients 
in  (EQ  34)  was  tested  and  verified.  In  addition,  the  dependence  of  the  material  properties 
on  temperature  will  be  discussed. 


u  rj  afliaffjHi  ms 


This  phase  in  any  program  intensive  application  is  critical  because  only  through  care¬ 
ful  checking  and  the  use  of  multiple  test  procedures  can  one  be  sure  that  the  program  is 
doing  what  it  was  designed  to  do.  Typically,  a  known  problem  is  solved  by  some  method 
other  than  the  method  used  in  the  program  and  compared  to  the  results  of  the  program.  In 
this  section,  verifications  with  hand  calculations  and  verification  with  an  exact  integration 
program  will  be  presented. 

3.1.1.  Verification  With  Hand  Calculations 

The  program  “Polystress”  was  checked  by  hand  for  a  second  order  polynomial 
approximation  of  the  integrand  in  (EQ  34).  For  a  second  order  polynomial  the  [B]  matrix 
is: 

00200  2£  0  0  2£2 

[fi]  =  0  0  000  0  2  2r\  2ti2  ^Q35) 

0  0  0  0  1  2ti  0  2C  4Cn 

Note  that  [B]  is  a  3xN  where  (N=(n+1)2).  The  transpose  of  [B]  was  multiplied  by  an 
assumed  inverted  extensional  stiffness  matrix  [a].  The  resulting  matrix  was  then  multi¬ 
plied  with  [B]  and  exactly  integrated  over  the  interval  -1  to  1.  This  matrix  was  then  com¬ 
pared  to  the  results  of  numerically  integrating  the  same  equation  and  they  were  identical  to 
the  sixth  significant  figure.  The  same  procedure  was  followed  for  the  thermal  load  vector 
on  the  right  hand  side  of  (EQ  34)  and  again  the  results  were  identical  to  the  sixth  signifi¬ 
cant  figure. 

3.1.2  Comparison  with  Exact  Integration  Program 

This  program  was  written  as  a  supplement  for  an  exact  integration  code  that  would 
calculate  not  only  the  in-plane  resultant  forces  but  also  the  out-of-plane  resultant  moments 
for  a  thermally  loaded  fibrous  composite  plate.  At  the  present  time  the  code  was  exactly 
integrating  the  extensional  stiffness  matrix  with  no  temperature  dependent  relations 
included.  Thus,  I  could  compare  my  results  with  the  exact  integration  program.  Clearly, 
the  hand  calculations  and  the  results  from  the  numerical  integration  could  be  used  to  ver¬ 
ify  the  exact  integration  program  as  well.  As  it  turned  out  all  three  methods  produced  the 
same  results  and  thus  the  numerical  and  exact  integration  phases  of  the  respective  pro¬ 
grams  had  been  verified. 


3.2  TEMPERATURE  DEPENDENT  MATERIAL  PROPERTIES 

The  modulus  of  elasticity,  E,  and  the  coefficient  of  thermal  expansion,  a,  are  not  tem¬ 
perature  independent.  In  data  given  by  (REF  8)  the  modulus  of  elasticity  decreases  almost 
linearly  with  increasing  temperature  and  the  coefficient  of  thermal  expansion  increases 
almost  linearly  with  increasing  temperature.  References  9  and  10  also  use  this  type  of  lin¬ 
ear  relationship  for  their  temperature  dependent  properties. 

In  (REF  10)  the  temperature  dependence  of  the  modulus  of  elasticity  is  expressed  in 
the  form 

E,{T)  =  £?(l-yT)  =  £j(l-f5(£)) 

1 1 

E2(T)  =  E°2(\-YT)  =  E2(  1-(H^)) 

1  (EQ  36) 

where  p=YT i  and  E°i  and  E°2  are  the  values  of  the  moduli  of  elasticity  along  the  1  and  2 
directions  at  the  reference  temperature  T0.  Here  T  denotes  the  temperature  excess  above 
the  reference  temperature  at  any  point  and  (J  denotes  the  percent  change  of  the  modulus  of 
elasticity  over  the  specified  temperature  interval  (T0<T<T  j)  Values  of  y  were  calculated 
for  a  25%,  and  50%  reduction  in  the  modulus  of  elasticity.  These  same  values  were  used 
for  the  corresponding  increase  of  the  thermal  expansion  coefficient.  For  changes  of  25% 
and  50%  the  values  for  y  are:  0.0005,  and  0.001  respectively,  forT0=75  and  Ti=550 
degrees  Fahrenheit. 

Since  the  variational  method  allows  for  a  continuous  solution,  a  quadratic  curve  fit 
which  very  closely  matches  an  experimental  temperature  distribution  obtained  from 
NASA-Dryden  was  employed  (REF  11).  Figure  1  displays  the  temperature  distribution 
used  for  this  study  which  is  given  by  the  relation  T(x)=  dT(x)  +  75.0.  A  one  dimensional 
curve  fit  was  accomplished  due  to  time  constraints  and  because  this  is  a  parametric  study 
to  give  baseline  theoretical  results  on  very  specific  problem  to  be  used  for  a  more  general 
problem  at  NASA.  The  curve  fit  used  in  the  program  is  from  Reference  12. 
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FIGURE  1.  Computational  Temperature  Distribution 


The  materials  scheduled  for  testing  were  Graphite/Epoxy,  Boron/Epoxy,  E-Glass/ 
Epoxy,  and  Aluminum.  Properties  for  these  materials  were  taken  from  References  2  and 
8.  These  materials  have  certain  temperature  limitations.  However,  since  this  was  a  pre¬ 
liminary  parametric  study  of  the  trends  for  a  thermally  loaded  composite  plate  this  limita¬ 
tion  was  overlooked.  The  materials  selected  could  have  been  changed,  or  the  temperature 
distribution  reduced,  but  the  materials  are  typical  for  test  cases  studied  in  most  applica¬ 
tions. 

3.3  CONVERGENCE 

In  most  approximation  routines  it  is  necessary  to  show  that  the  computer  code  is  con¬ 
verging  to  a  single  solution,  whether  or  not  this  is  a  true  solution  depends  on  the  method 
and  comparison  with  experimental  data.  For  the  present  study,  the  solution  should  con¬ 
verge  as  the  order  of  the  polynomial  approximation  increases. 

For  the  present  study,  convergence  was  verified  by  taking  the  square  root  of  the  sum  of 
squares  for  the  in-plane  loads  at  a  specified  number  of  points  within  the  plate  and  sum¬ 
ming  them  for  the  entire  plate.  This  process  was  repeated  for  different  order  polynomials. 
Since  this  method  is  dependent  on  the  number  of  points,  the  same  number  of  points  was 
used  for  all  the  test  cases:  77  grid  points,  with  11  in  the  x-direction  and  7  y  locations  for 
each  x  location .  Contour  plots  for  eighth  and  tenth  order  approximating  polynomials  are 
shown  in  Attachment  1.  Clearly,  the  results  for  the  two  different  order  of  polynomials  are 
very  close  and  would  appear  to  support  the  premise  that  the  solution  is  converging.  How¬ 
ever,  additional  orders  of  polynomial  approximations  were  done  to  verify  that  conver¬ 
gence  had  been  reached.  Using  the  sum  of  the  square  root  of  the  sum  of  squares,  a 
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percent  difference  of  2.5  was  calculated  for  the  entire  plate  for  polynomial  approximations 
of  8  and  10  respectively.  The  difference  between  an  eighth  and  tenth  order  approximation 
at  the  center  of  the  plate,  x=25  inches  and  y=6  inches,  is  approximately  2%. 


Additional  convergence  test  cases  were  accomplished  for  higher  order  polynomials. 
However,  the  coefficient  matrix  [A]  which  is  equivalent  to  [B]T[a][B]  in  (EQ  35) 
becomes  ill-conditioned  upon  inversion,  thus  the  accuracy  is  reduced.  The  difference 
between  a  tenth  and  twelfth  order  approximation  was  approximately  1.8%,  and  the  differ¬ 
ence  between  higher  order  polynomials  was  slightly  worse  because  the  [A]  matrix  became 
too  ill-conditioned.  A  tenth  order  polynomial  approximation  was  chosen  for  the  test  cases 
because  it  was  the  best  trade-off  between  time  and  percent  convergence. 


The  results  from  the  contour  plots  closely  matched  experimental  data  from  Reference 
11  on  similar  test  cases.  This  analysis  has  shown  the  verification  of  the  program 
“Polystress”.  Finally,  the  program  has  solved  a  problem  with  a  known  solution  quite 
accurately,  and  it  can  now  be  applied  to  problems  whose  solutions  are  unknown  with  some 
validity. 

3.4  TEST  MATRIX 


Typically  a  test  matrix  is  set  up  prior  to  the  test  period  for  any  project.  This  is  to 
ensure  that  the  actual  testing  time  is  used  efficiently.  If  something  peculiar  is  found  while 
processing  the  data  additional  testing  can  be  accomplished,  but  usually  only  after  the  ini¬ 
tial  matrix  is  completed. 

The  test  matrix  for  this  problem  is  shown  in  Table  1.  The  tests  in  Table  1  will  be  per¬ 
formed  for  a  fibrous,  laminated,  composite  plate  with  dimensions:  length=  50  in,  width= 
12  in,  thickness=.19  in.  A  tenth  order  polynomial  will  be  used  to  approximate  the  solu¬ 
tion,  and  an  eighth  order  polynomial  will  be  used  to  verify  convergence  and  input  as 
required. 
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TABLE  1.  Test  Matrix 


Material 

LavuD 

#  Lavers 

Gamma 

Theta  Ranee 

Graphite/Epoxy 

[0] 

1 

0.00 

.... 

Boron/Epoxy 

[0] 

1 

0.00 

.... 

E-Glass/Epoxy 

[0] 

1 

0.00 

— 

Aluminum 

[03 

1 

0.00 

— 

Graphite/Epoxy 

[61 

1 

0.0005 

0  to  90 

Boron/Epoxy 

[6] 

1 

0.0005 

0  to  90 

E-Glass/Epoxy 

[61 

1 

0.0005 

0  to  90 

Aluminum 

[0] 

1 

0.0005 

— • 

Graphite/Epoxy 

[+0/-0/+O] 

3 

0.0005 

0  to  90 

Boron/Epoxy 

[+e/-e/+8] 

3 

0.0005 

0  to  90 

E-Glass/Epoxy 

[+0/-0/+81 

3 

0.0005 

0  to  90 

E-Glass/Epoxy 

[+0/-0/+0/-j 0/+0]  5 

0.0005 

0  to  90 

Graphite/Epoxy 

[0] 

1 

0.001 

— - 

Boron/Epoxy 

[03 

1 

0.001 

.... 

E-Glass/Epoxy 

[01 

1 

0.001 

.... 

Aluminum 

[0] 

1 

0.001 

— - 

Here  y  is  the  factor  applied  to  the  modulus  of  elasticity  and  the  coefficient  of  thermal 
expansion.  The  cases  investigated  are  all  for  symmetric  laminates  since  there  is  no  cou¬ 
pling  between  extensional  and  bending  stiffnesses  for  this  layup. 

Two  types  of  symmetric  laminates  are  studied:  1)  laminates  with  (single/multiple)  spe¬ 
cially  orthotropic  layers,  and  2)  laminates  with  multiple  generally  orthotropic  layers. 
Specially  orthotropic  implies  that  the  stiffness  matrix  [Q]  and  the  reduced  stiffness  matrix 
[Q]  are  identical,  (i.e.  the  principal  material  coordinates  of  the  lamina  is  aligned  with  glo¬ 
bal  laminate  axis  for  the  plate)  (REF  2).  Item  2  contains  the  special  classification:  regular 
symmetric  angle-ply  laminate.  This  specification  refers  to  orthotropic  lamina,  of  equal 
thicknesses,  with  opposite  signs  of  the  angle  of  orientation  of  the  principal  material  axis 
relative  to  the  laminate  axis  system  for  the  problem. 
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a 

The  primary  purpose  of  this  study  is  to  perform  a  parametric  investigation  of  the  prin¬ 
cipal  in-plane  force  resultants  for  a  symmetrically  laminated,  50”x  1 2”x0. 1 9”coi  nposite, 
plate  for  a  thermal  loading  condition.  The  boundary  conditions  for  the  middle  surface 
stresses  for  this  study  are  free-free.  This  condition  implies  that  both  the  midsurface  stress 
components  normal  to  the  boundary  and  tangential  to  the  boundary  are  zero.  This  is  the 
conventional  stress  configuration  (REF  6).  Four  different  materials  were  investigated. 


TABLE  2.  Material  Properties 


Property 

Graphite/Epoxy 

Boron/Epoxy 

E-Glass/Epoxy 

Aluminum 

Ei 

30e6  psi 

30c6  psi 

7.8c6  psi 

10.3c6  psi 

e2 

0.75c6  psi 

3.0c6  psi 

2.6c6  psi 

10.3c6  psi 

v12 

0.25 

0.30 

0.25 

0.33 

Gi2 

1.3c6  psi 

1.0e6  psi 

,375c6  psi 

Isotropic 

<*1 

-0.21e-6/F° 

3.5e-6  /F° 

3.5e-6/F° 

12.8e-6/F° 

<*2 

16.0c-6  /F° 

11.4c  6 /F° 

11.4c-6/F° 

12.8c-6/F° 

Properties  were  taken  from  References  1,2,  and  8. 


4.1  COMPARISON  OF  PRINCIPAL  IN-PLANE  FORCE  RESULTANTS 

Variations  in  lamina  orientation  angles  and  number  of  layers  are  considered  for  the 
free-free,  laminated,  composite,  plate.  The  principal,  in-plane  force  resultants  for  a  single¬ 
lamina  and  a  three-layered,  regularly  symmetric  angle  ply  laminate  as  a  function  of  orien¬ 
tation  angle  are  shown  in  Figures  2,  and  3  respectively.  Principal  in-plane  force  resultant, 
is  essentially  the  maximum  in-plane  force  resultant  for  the  plate  regardless  of  whether  it  is 
Nx,  Ny,  Nxy.  In  addition,  the  positive  and  negative  principal  in-plane  force  resultants 
were  graphed  separately  to  aid  in  analysis. 
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FIGURE  3.  Principal  (+)  in-plane  force  resultant  for  a  3-layered  laminate. 


In  both  Figures  2  and  3  the  magnitude  of  the  principal  resultant  force  is  greatest  at  0,  and 
conversely  at  180  degrees.  This  is  a  result  of  the  one  dimensional  thermal  loading  applied 
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to  the  plate.  When  the  lamina  is  aligned  with  the  load  the  in-plane  stresses  are  the  greatest 
because  the  modulus  of  elasticity  for  the  lamina  is  aligned  with  the  load  instead  of  being 
at  an  angle  relative  to  the  load. 

Furthermore,  the  relation  between  the  material  properties,  Table  1 ,  and  the  relations 
expressed  in  Figure  2  is  pronounced.  Initially,  Boron/Epoxy  has  a  higher  force  resultant 
than  any  of  the  other  materials  because  the  moduli  of  elasticity  for  Boron/Epoxy  are  as 
high  or  higher  than  any  of  the  other  materials.  Also  as  the  orientation  angle  exceeds 
approximately  45  degrees  the  force  resultants  for  Boron/Epoxy  and  E-Glass/Epoxy  are 
nearly  identical.  This  result  is  due  to  the  same  coefficients  of  thermal  expansion  being 
used  for  both  materials,  and  the  relative  decline  in  importance  of  Ej  and  the  increasing 
importance  of  E2  as  the  orientation  angle  approaches  90  and  E2  is  aligned  in  the  x-direc- 
tion  of  the  laminate.  Note  the  similarity  in  the  values  for  E2  for  both  materials.  The  iso¬ 
tropic  case  for  aluminum  was  shown  on  both  graphs  as  a  reference  for  the  other  values. 
Tabular  values  for  this  case  and  the  other  test  cases  are  in  Attachment  2.  Similar  results  to 
Figures  2  and  3  are  given  in  Reference  2. 

Figure  3  illustrates  the  same  trends  as  Figure  2,  however,  the  magnitude  of  the  princi¬ 
pal  resultant  forces  is  greater  at  ten  degrees  than  at  zero  for  all  the  composite  materials 
except  E-Glass/Epoxy.  Another  interesting  observation  involves  the  comparison  of  the 
unidirectional  lamina  at  different  orientations  in  Figure  2  with  the  angle-ply  data  in  Figure 
3.  For  angles  larger  than  0°  but  less  than  45°,  the  angle-ply  has  about  a  50  to  60  percent 
higher  resultant  force  than  does  the  unidirectional  lamina.  However,  at  approximately  40 
to  45  degrees  and  higher,  the  single  lamina  has  a  greater  resultant  force.  Again,  E-Glass/ 
Epoxy  exhibits  the  same  trends  but  on  an  order  of  magnitude  lower.  The  trends  for  higher 
loads  at  10°  than  at  0°,  and  the  differences  between  the  data  for  the  single  ply  and  the  data 
for  the  multi-layer  angle  ply,  are  the  results  of  mechanical  and  thermal  interactions 
between  the  layers  (REF  2). 

The  results  for  Graphite/Epoxy  are  more  dramatic  than  any  of  the  other  cases.  This  is 
because  the  coefficient  of  thermal  expansion  in  the  1  principal  direction  is  less  than  zero 
(REF  1).  This  implies  that  the  expansion  in  the  2  principal  direction  is  large  enough  to 
cause  a  contraction  in  the  1 -direction.  This  reaction  produces  thermal  reactions  between 
the  layers  which  causes  the  resultant  forces  to  be  larger  than  normal.  As  the  angle  is 
increased,  the  contribution  of  the  (*2  term  is  reduced,  and  the  relatively  low  value  of  elas¬ 
tic  moduli  in  the  2-direction  causes  the  rapid  decline  of  the  resultant  forces. 

The  relative  invariant  behavior  of  E-Glass/Epoxy  is  a  result  of  the  relatively  small 
change  between  the  elastic  moduli  in  the  1  and  2  principal  material  directions.  Compared 
to  the  other  materials,  except  of  course  Aluminum,  the  difference  between  Ej,  and  E2,  is 
order  of  magnitudes  lower.  This  difference  explains  the  results  in  Figures  2  and  3. 

Figures  4  and  5  illustrate  the  principal  negative  force  resultants  for  the  same  layup  as 
Figures  2  and  3. 


0.00 


FIGURE  4.  Principal  (•)  in-plane  force  resultant  for  a  single  lamina. 


The  trends  for  the  data  in  Figure  4  are  the  same  as  those  presented  for  Figure  1. 


FIGURE  5.  Principal  (-)  in-plane  force  resultant  for  a  3-layered  laminate. 


The  data  in  Figures  4  and  5  follows  the  same  trends  as  described  earlier  for  the  principal, 
positive,  force  resultant. 


A  comparison  of  the  effect  of  the  number  of  layers  on  the  principal  in-plane  force 
resultants  was  accomplished  for  the  free-free,  E-Glass/Epoxy,  laminated  plate.  Figures  6, 
and  7  show  the  results  for  1, 3,  and  5  layer  angle  ply  laminates  as  a  function  of  lamina  ori¬ 
entation  angle. 

As  discussed  earlier  for  the  comparison  between  the  single  lamina  in  Figure  2  and  the 
3-layered  laminate  in  Figure  3,  the  interlaminar  mechanical  and  thermal  interactions 
account  for  the  multi-layered  laminate  having  larger  in-plane  resultant  forces  than  the  uni¬ 
directional  lamina  for  orientation  angles  greater  than  0°  and  less  than  approximately  45° 
degrees.  For  angles  greater  than  45°  the  single-layered  lamina  has  higher  force  resultants. 
The  interlaminar  effects  explanation  is  also  valid  for  Figures  6  and  7.  From  0°  to  approxi¬ 
mately  40°  the  resultant  forces  for  the  3  and  5  layer  laminates  are  higher  than  those  for  the 
single  layer  case.  For  angles  greater  than  40°,  the  single-layer  forces  are  greater.  Also  the 
difference  between  the  force  resultants  for  the  3  and  5  layer  laminates  is  much  less  than 
that  for  the  1  and  3  layer  laminates.  One  would  expect  the  difference  for  each  successive 
layer  to  get  smaller  and  smaller  based  on  these  assumptions. 

To  verify  this  assumption  a  test  case  was  ran  for  a  [20/-20/20/-20/20/-20/20]  E-Glass/ 
Epoxy  laminate.  The  principal  in-plane  force  resultants  are:  45.129  IbAn  and  -92.251  lb/in 
respectively.  The  previous  results  for  the  E-Glass/Epoxy  1 -layer,  3-layer,  and  5-layer 
laminate  are:  39.998  &  -85.159,  44.433  &  -91.055,  and  45.030  &  -92.128  lb/in  respec¬ 
tively.  These  results  are  shown  graphically  in  Figure  8. 


n=10 

thickness=.19  in 
length=  50  in 
width=  12  in 
¥=.0005 

— 0 —  1  Layer 

— □ —  3  Layer 

\o  — & —  5  Layer 


Theta  (degrees) 


FIGURE  6.  E-Glass/Epoxy  multiple  layer  comparison  (+). 
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FIGURE  8.  Effect  of  multiple  layers  on  the  force  resultants  for  E-Glass/Epoxy. 


I 
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Recall  that  y  is  the  factor  applied  to  the  moduli  of  elasticity  and  the  coefficient  of  ther¬ 
mal  expansion  to  account  for  temperature  dependency.  Figure  9  displays  the  results  for 
the  investigation  of  the  four  different  materials  as  specially  orthotropic  single  lamina. 

The  magnitudes  of  the  principal  in-plane  force  resultants  increased  as  y  increased. 
However,  according  to  Hooke’s  Law,  (EQ  1),  the  stress  is  directly  related  to  the  moduli  of 
elasticity.  Thus  it  would  seem  that  the  stress  should  decrease  with  an  increase  in  y  since  y 
is  the  temperature  degradation  factor  applied  to  the  moduli  of  elasticity.  But,  y  is  also 
applied  to  the  coefficient  of  thermal  expansion,  a,  which  increases  with  temperature. 
Therefore,  the  effects  of  a  on  the  thermal  load  vector  must  be  great  enough  to  overcome 
the  effects  of  y  on  the  moduli  of  elasticity. 


— A —  Graphite/Epoxy 
— B —  Boron/Epoxy 
— © —  E-Glass/Epoxy 
180.00.  Aluminum 

16°.°°  7  JPr^^ 

140.00  r 


§  100.00, 

K 

£  80.00 

3-  60.00 


0.00025  0.0005  0.00075  0.001 

Gamma 


_  FIGURE  9.  Effects  of  gamma  on  principal  in-plane  force  resultant. 


4.4  CONTOUR  PLOTS  FOR  THE  FOUR  DIFFERENT  MATERIALS 


Contour  plots  for  the  four  different  materials  investigated  during  this  study  are  pre¬ 
sented  in  Attachment  3.  The  plots  shown  are  the  force  resultant  in  the  x-direction,  Nx,  for 
zero  y,  and  are  presented  to  give  the  basic  characteristics  of  each  material  in  a  graphical 
instead  of  tabular  method.  In  addition,  the  contour  plot  for  each  material  is  symmetrical 
because  the  gradient  of  the  temperature  distribution  used  for  this  study  was  symmetrical. 

The  contour  plot  for  Graphite/Epoxy  exhibits  the  largest  range  of  force  resultants  for 
any  of  the  materials,  and  E-Glass/Epoxy  exhibits  the  smallest  range.  These  results  are 
functions  of  the  material  properties  given  in  Table  2. 
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An  analytical  approach  for  determining  the  in-plane  force  resultants  for  symmetrically 
laminated,  composite  plates,  subjected  to  a  one  dimensional  thermal  loading  has  been  pre¬ 
sented.  The  Galerkin  method  was  used  to  approximately  solve  the  compatibility  equa¬ 
tion.  A  linear  response  of  the  material  properties  with  temperature  was  employed.  The 
results  for  the  baseline  case  of  aluminum  for  this  study  were  within  1  percent  of  an  exact 
solution  for  the  same  case.  In  addition,  the  results  demonstrate  the  dependency  of  the  in¬ 
plane  force  resultants  with  orientation  angle,  and  the  relation  between  single-layer  laminas 
and  multi-layer  laminates.  The  inter-lamina  effects  for  the  multi-layered  laminate  resulted 
in  higher  force  resultants  up  to  an  orientation  angle  of  45  degrees  as  compared  to  the  force 
resultants  for  the  unidirectional  laminate.  The  temperature  dependency  of  the  engineer¬ 
ing  constants  resulted  in  higher  force  resultants  for  an  increase  in  y,  the  temperature 
dependency  factor.  These  investigations  should  provide  an  excellent  database  for  experi¬ 
mental  studies  of  composite  plates  at  NASA-Dryden.  Further  studies  will  investigate  out- 
of-plane  motion,  for  a  two  dimensional  temperature  distribution  using  exact  integration 
techniques. 
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Attachment  2.  Table  of  Results 


bl 


Principal  In-Plane  Force  Resultant 

Length=50  in,  Width=12  in.  Thickness^  .19  in,  n= 10 


Case# 

Material 

Layup 

Gamma 

Pos 

Neg 

1 

Graphite/Epoxy 

[0] 

0.00 

95.078 

-297.794 

2 

Graphite/Epoxy 

10] 

0.0005 

123.255 

-389.614 

3 

Graphite/Epoxy 

[0] 

0.001 

141.107 

-458.437 

4 

Boron/Epoxy 

[0] 

0.00 

116.701 

-264.507 

5 

Boron/Epoxy 

[0] 

0.0005 

152.819 

-365.965 

6 

Boron/Epoxy 

[0] 

0.001 

177.907 

-434.965 

7 

E-Glass/Epoxy 

[0] 

0.00 

40.120 

-80.860 

8 

E-Glass/Epoxy 

[0] 

0.0005 

56.615 

-118.856 

9 

E-Glass/Epoxy 

[0] 

0.001 

67.984 

-143.959 

10 

Aluminum 

[0] 

0.00 

60.251 

-121.167 

11 

Aluminum 

[0] 

0.0005 

92.166 

-191.157 

12 

Aluminum 

[0] 

0.001 

111.402 

-229.105 

13 

Graphite/Epoxy 

[10] 

0.0005 

93.460 

-185.353 

14 

Graphite/Epoxy 

[20] 

0.0005 

37.331 

-71.444 

15 

Graphite/Epoxy 

[30] 

0.0005 

15.573 

-32.953 

16 

Graphite/Epoxy 

[40] 

0.0005 

9.033 

-17.200 

17 

Graphite/Epoxy 

[50] 

0.0005 

3.989 

-9.216 

18 

Graphite/Epoxy 

[60] 

0.0005 

2.366 

-4.844 

19 

Graphite/Epoxy 

[70] 

0.0005 

0.936 

-2.327 

20 

Graphite/Epoxy 

[80] 

0.0005 

0.147 

-0.432 

21 

Graphite/Epoxy 

[90] 

0.0005 

0.340 

-0.109 

22 

Boron/Epoxy 

[10] 

0.0005 

106.450 

-245.463 

23 

Boron/Epoxy 

[20] 

0.0005 

61.323 

-120.464 

24 

Boron/Epoxy 

[30] 

0.0005 

34.410 

-66.828 

b2 


Case# 

Material 

Layup 

Gamma 

Pos 

Neg 

25 

Boron/Epoxy 

[40] 

0.0005 

27.406 

41.508 

26 

Boron/Epoxy 

[50] 

0.0005 

15.297 

-30.075 

27 

Boron/Epoxy 

[60] 

0.0005 

11.202 

-23.153 

28 

Boron/Epoxy 

[70] 

0.0005 

8.803 

-19.325 

29 

Boron/Epoxy 

[80] 

0.0005 

7.520 

-18.562 

30 

Boron/Epoxy 

[90] 

0.0005 

7.197 

-18.305 

31 

E-Glass/Epoxy 

[10] 

0.0005 

51.472 

-106.982 

3? 

E-Glass/Epoxy 

[20] 

0.0005 

39.998 

-85.159 

33 

E-Glass/Epoxy 

[30] 

0.0005 

28.629 

-59.117 

34 

E-Glass/Epoxy 

[40] 

0.0005 

21.013 

41.738 

35 

E-Glass/Epoxy 

[50] 

0.0005 

15.249 

-30.594 

36 

E-Glass/Epoxy 

[60] 

0.0005 

11.183 

-22.900 

37 

E-Glass/Epoxy 

[70] 

0.0005 

8.500 

-17.514 

38 

E-Glass/Epoxy 

[80] 

0.0005 

6.959 

-15.043 

39 

E-Glass/Epoxy 

[90] 

0.0005 

6.455 

-14.364 

40 

Graphite/Epoxy 

[10/-10/10] 

0.0005 

189.613 

433.888 

41 

Graphite/Epoxy 

[20/-20/20] 

0.0005 

134.863 

-263.107 

42 

Graphitc/Epoxy 

[30/-30/30] 

0.0005 

41.275 

-78.245 

43 

Graphitc/Epoxy 

[40/40/40] 

0.0005 

6.474 

-11.168 

44 

Graphitc/Epoxy 

[50/-50/50] 

0.0005 

3.513 

-1.765 

45 

Graphitc/Epoxy 

[60/-60/60] 

0.0005 

7.723 

-2.227 

46 

Graphitc/Epoxy 

[70/-70/70] 

0.0005 

4.835 

-1.266 

47 

Graphitc/Epoxy 

[80/-80/80] 

0.0005 

1.386 

-0.378 

48 

Graphite/Epoxy 

[90/-90/90] 

0.0005 

0.340 

-0.109 

49 

Boron/Epoxy 

[10/-10/10] 

0.0005 

165.984 

-365.550 

50  Boron/Epoxy  [20/-20/20]  0.0005  135.700  -279.620 


Boron/Epoxy  [30/-30/30] 


63 


51 


0.0005 


71.374 


-141.720 


Case  # 

Material 

Layup 

Gamma 

Pos 

Neg 

52 

Boron/Epoxy 

[40/  -40/40] 

0.0005 

26.095 

47.760 

53 

Boron/Epoxy 

[50/-50/50] 

0.0005 

7.792 

-14.778 

54 

Boron/Epoxy 

[60/-60/60] 

0.0005 

4.753 

-13.693 

55 

Boron/Epoxy 

[70/-70/70] 

0.0005 

5.826 

-18.869 

56 

Boron/Epoxy 

[80/-80/80] 

0.0005 

6.962 

-20.227 

57 

Boron/Epoxy 

[90/-90/90] 

0.0005 

7.197 

-18.304 

58 

E-Glass/Epoxy 

[10/-10/10] 

0.0005 

53.408 

-111.120 

59 

E-Glass/Epox> 

[20/-20/20] 

0.0005 

44.433 

-91.055 

60 

E-Glass/Epoxy 

[30/-30/30] 

0.0005 

31.244 

o4.359 

61 

E-Glass/Epoxy 

[40/40/40] 

0.0005 

19.760 

39.684 

62 

E-Glass/Epoxy 

[50/-50/50] 

0.0005 

11.991 

-24.000 

63 

E-Glass/Epoxy 

[60/-60/60] 

0.0005 

8.193 

-16.454 

64 

E-Glass/Epoxy 

[70/-70/70] 

0.0005 

6.853 

-14.559 

65 

E-Glass/Epoxy 

[80/-8G/80] 

0.0005 

6.494 

-14.346 

66 

E-Glass/Epoxy 

[90/-90/90] 

0.0005 

6.421 

-14.289 

67 

E-Glass/Epoxy 

[10/-10/10/-10/10] 

0.0005 

53.866 

-112.070 

68 

E-Glass/Epoxy 

[20/-20/20/-20/20] 

0.0005 

45.030 

-92.128 

69 

E-Glass/Epoxy 

[30/-30/30/-30/30] 

0.0005 

31.627 

-64.973 

70 

E-Glass/Epoxy 

[40/40/40/40/40] 

0.0005 

19.790 

-39.812 

71 

E-Glass/Epoxy 

[50/-50/50/-50/50] 

0.0005 

11.849 

-23.661 

72 

E-Glass/Epoxy 

[60/-60/60/-60/60] 

0.0005 

8.047 

-16.102 

73 

E-Glass/Epoxy 

[70/-70/70/-70/70] 

0.0005 

6.770 

-14.471 

74 

E-Glass/Epoxy 

[80/-80/80/-80/80] 

0.0005 

6.495 

-14.375 

75 

E-Glass/Epoxy 

[90/-90/90/ -90/90] 

0.0005 

6.455 

-14.364 

76 

E-Glass/Epoxy 

[25/-25/251 

0.0005 

38.807 

-78.003 

77 

E-Glass/Epoxy 

[32/32/32] 

0.0005 

28.533 

-58.944 

78 

E-Glass/Epoxy 

[33/-33/33J 

0.0005 

27.333 

-56.304 

b4 


Case  # Material_ Layup_ Gamma_ Pos_ Neg 

79  _ E-Glass/Epoxy  [45/45/45] _ 0.0005 _ 15.340  -30.713 

80  _ E-Glass/Epoxy  [65/-65/6S] _ 0.0005 _ 7.325 _ -14,921 

81  _ E-Glass/Epoxy  [75/-75/75J _ 0.0005 _ 6.612 _ -14.413 

82  _ E-Glass/Epoxy  [957-95/95] _ 0.0005 _ 6.438 _ -14,305 

83  _ E-Glass/Epoxy  [110/-110/110]  0.0005 _ 6.853 _ -14,559 

84  _ E-Glass/Epoxy  [125/-125/125]  0.0005 _ 9£71 _ -19.382 

85  _ E-Glass/Epoxy  [135/-135/135]  0.0005 _ 15.340 _ 30.713 

86  _ E-Glass/Epoxy  [145/- 145/145}  0,0005 _ 25.059 _ -51.184 

87  _ E-Glass/Epoxy  [150/- 150/1 50]  0.0005 _ 31.244 _ -64,339 

88  _ E-Glass/Epoxy  [165/-165/1651  0-0005 _ 49.659 _ -102.502 

89  _ E-Glass/Epoxy  [175/- 175/1 75]  0.0005  55.603  -116.328 

E-Glass/Epoxy  [180/- 180/1 80] 


90 


0.0005 


56.615 


-118.856 


Attachment  3.  Contour  Plots  For  Different  Materials 
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Attachment  4.  Program  Listing 


PROGRAM  POLYSTRESS 

Q  ************************************************* 

C  **  POLYNOMIAL  RITZ  SOLUTION  FOR  IN  PLANE  STRESS  * 

C  **  EVALUATION  WITH  A  THERMAL  LOADING  CONDITION  * 
c  **  * 

C  **  THE  IN  PLANE  FORCE  RESULTANTS  ARE  WRITTEN  TO  * 

C  **  FILE  CALLED  SDISTRIBUTION.DAT  * 

c  **  * 

c  **  Darren  Knipp  * 

0  ***★*****■*******★***★*★*★*★*★★****★***★★★*★★★★*** 
implicit  real*8  (a-h,o-z) 
parameter (in=200) 

dimension  am{in,in),  bm(in,in),  cm(in,in),  av(in),  bv(in) 

dimension  bmat (in, in) ,bmgrid(in, in) 

dimension  amd (in, in) , avd (in) , tk (25) , theta (25) 

dimension  x(9),T(9),TE(9), break (9) , tcoef (4,9) 

common  /worksp/rwksp 

real  rwksp (2*in**2+3*in) 

call  iwkin (2*in**2+3*in) 

c  **  open  output  data  files  *********** 

open (unit=9, file-' extstif f . dat' , status-' unknown' ) 
open (unit-10, file=' thermload.dat' ,  status*' unknown' ) 
open (unit=ll, file*' stress f . dat' , status*' unknown' ) 
open (unit=12, file*' distrib.dat', status*'  unknown'  ) 

c  ******  Temperature  data  *********** 
c  ** 

C  **  dT (x)  =  0.2*x~2,  T(x)=dT(x)+75. *********** 

data  x/0., 6. 25, 12. 5, 18. 75, 25., 31. 25, 37. 5, 43. 75, 50./ 

data  T/75 . , 82.81,106.25, 145 . 31, 200 ., 270 . 31 , 356 . 25, 457.81,575./ 

data  TE/75., 82. 81, 106. 25, 145. 31, 200. ,270.31,356.25,457.81, 575./ 

C  ***********  DATA  INPUT  ************** 

call  input (a, b, h, ip, nf , nx, ny, nl, tk, theta, t factor, gfac) 

c  **  scale  the  room  temperature  wrsp  to  room  temp,  of  75  deg.  ** 
do  5  i=l, 9 

5  T (i) =t fact or* (T (i) -75 . ) +75. 

c  **  curve  fit  the  temperature  data  ** 
call  dcsakm(9, x,T, break, tcoef ) 

c  ***********  COMPUTE  THE  STRESS  FUNCTION  COEFFICIENTS  *** 

call  sfunction (a,b,h,nf,tk, theta, nl, am, bm, cm, av, bv, bmat, 

+  amd, avd, break, tcoef , ip, gfac) 

c  **  calculate  the  stress  distribution  and  write  to  output  file 

call  sdistribution (a, b, nf , nx,  ny,  av,  bmgrid) 

stop 

end 

c  ************************************************************ 
subroutine  input (a, b, h, ip, nf , nx, ny, nlayer, tlayer, theta, 

+  tf actor, gf actor) 

c  ** 

c  **  This  subroutine  prompts  the  user  for  the  problem  information, 
c  ** 

c  ************************************************************** 
implicit  real*8 (a-h, o-z) 
character*30  namel,mat 
dimension  tlayer (25) , theta (25) 
pi-3.141592653589793 

c  **  Define  the  problem  ** 


print*,'  Define  material  properties:  1,  Graphite/Epoxy 


d2 


2,  Boron/Epoxy' 

3,  E-Glass/Epoxy 

4,  Aluminum' 


print*, ' 
print*, ' 
print*, ' 
read (5, *) ip 
print*,'  ' 

print*,'  Input  the  output  data  filename  within  single  quotes  ' 
read (5,*)  namel 

open (unit=3, file=namel, status='  unknown' ) 
print*,'  ' 

print*, '  Input  order  of  polynomial  for  stress  function' 
read(5,*)  nf 

print*,'  Input  the  plate  length  and  width  ' 
read(5,*)  a,b 

print*, ' Input  the  number  of  layers:  nl  (nl  max=25)  ' 
read (5,*)  nlayer 

print*, 'For  each  layer  input:  thickness  and  theta  ' 
do  10  i=l, nlayer 
read(5,*)  tlayer(i) , theta (i) 
continue 

print*, ' Input  factor  to  apply  to  temperature  distribution' 
print*,' (0.0  =  room  temp.,  1.0  =  actual  lab  temp.)' 
read(5,*)  tfactor 

print*, ' Input  factor  to  apply  to  property  variance  with  temp.' 
print*,' (0<=  gfactor  =<1)' 
read(5,*)  gfactor 

print*, 'Dimension  of  your  grid  for  the  stress  distribution :x, y 
read (5,  *)  nx, ny 

if  (ip  .eq.  1)  then 
mat=' Graphite/Epoxy' 
endif 

if  (ip  .eq.  2)  then 
mat=' Boron/Epoxy' 
endif 

if  (ip  .eq.  3)  then 
mat=' E-Glass/Epoxy' 
endif 

if  (ip  .eq.  4)  then 
mat=' Aluminum' 
endif 

echo  input  data  ** 

write (3, 100) 
write (3, 150)  mat 
write (3, *) 


write (3, *) ' 
write(3,*)  ' 
write (3, *) 


Input  Data 


write (3, *) 
write (3, 200)  nlayer 
write (3, 300) 
do  20  i=l, nlayer 

write (3, 400)  i, t layer (i) , theta (i) 
continue 

compute  total  thickness  and  transform  ply  angles  to  radians 
h=0. 0 

do  30  i=l, nlayer 
theta (i) =theta (i) *pi/180 
h=h+tlayer (i) 
continue 
write(3,450)  h 

return 


format  (/'  Composite  Laminate  Stiffness'/) 
format('  Material:  ',a20/) 
format ('  Number  of  Layers  =  ',13/) 
format ('  Laminate  Geometry  '/,'  Layer 


theta' ) 


400 

450 


format (i4,4x,2f9.4) 

format (/' Laminate  Thickness  = 


f  9 . 4 ) 


end 

Q  ************************************************** 

subroutine  sfunction (a, b, h, n, tk, theta, nl, rm, rmb, tm, rv, rvb, bmat, 
+  rmd, rvd, break, tcoef, ip, gfac) 


c 

★  * 

c 

★  ★ 

This  subroutine  calculates  the  stress  function 

c 

*  * 

★  * 

coefficients: 

Fij. 

c 

c 

*  * 

variables : 

c 

*  * 

input: 

c 

★  ★ 

a 

—  plate  length  x  direction 

c 

*  * 

b 

—  plate  length  y  direction 

c 

★  * 

n 

—  order  of  stress  function  polynomial 

c 

*  * 

h 

—  plate  thickness 

c 

*  * 

output : 

c 

*  ★ 

rm 

—  [B] trans* [a] * [B]  order  N*N 

c 

*  ★ 

rv 

—  row  vector  containing  the  coefficient: 

Q  ★*★***★**★***★*★★*★**■*★***★★*★***★★*★***★★★**★***★★* 

implicit  real*8  (a-h,o-z) 
dimension  rm( (n+1) * (n+1) , (n+1) * (n+1) ) 
dimension  rmd( (n+1) * (n+1) , (n+1) * (n+1) ) 
dimension  rv ( (n+1 ) * (n+1 ) ) 
dimension  rvd ( (n+1 ) * (n+1 ) ) 
dimension  tm( (n+1) * (n+1)  ,  (n-3) * (n-3) ) 
dimension  rmb ( (n-3) * (n-3) , (n-3) * (n-3) ) 
dimension  rvb  ( (n-3)  *  (n-3)  ) 

dimension  bmat (3, (n+1) * (n+1) ) , amat (3,3) , arav (3) , tk (25) 
dimension  point (10) , weight (10) ,break(9) , tcoef (4, 9) , theta (25) 
nn= (n+1) * (n+1) 
m=3 

nnb= (n-3) * (n-3) 
rmfac=a*b/4 
rvfac=-a*b/4 
ng=10 

c  **  initialize  rm  matrix  to  zero  ****** 

do  10  i=l,nn 
rv  (i) =0 . 0 

do  10  j=l,nn 
rm(i, j) =0. 0 
10  continue 

c  ****  call  subroutine  for  numerical  integration  points  and  weights 
call  gquad (point, weight) 

c  **  evaluate  matrices  at  integration  point  ***** 

do  30  i=l,ng 
zeta=point (i) 
wz=weight (i) 

call  astiff(a,zeta,tk,theta,nl, break, tcoef, ip, gfac, amat, amv) 
j-0 

do  30  j=l,ng 
eta=point ( j) 
we=weight ( j) 

call  bmatrix (a,b, n, zeta, eta, bmat) 
call  btab (amat, bmat, rmd, m, nn) 
call  bta (1, nn,m, amv, bmat,  rvd) 
ii=0 

do  20  ii=l,nn 

rv (ii) =rv (ii) +rvd (ii) *wz*we 
j  j=0 

do  20  jj=l,nn 

rm(ii,  j  j)  =*rm(ii,  j  j)  +rmd(ii,  j  j)  *wz*we 
20  continue 

30  continue 


44 


40 


do  40  i=l,nn 
rv (i) =rvfac*rv (i) 
do  40  j=l,nn 
rm{i, j) =rmfac*rm(i,  j) 
continue 

c  **  write  rm {stiff ness  matrix)  &  rv(load  vector)  to  output  data  file 

write (10, 100)  (rv (i) , i=l, nn) 
do  50  j=l,nn 

write(9,100)  (rm(i, j) , i=l, nn) 

50  continue 

100  format (6el2 . 5) 

c  **  solve  matrix  equation  including  the  linear  transformation 

call  transform (n, tm) 

call  btab  (rm,  tm,  rmb,  nn,  nnb) 

call  bta (1, nnb, nn, rv,tm, rvb) 

call  dlslsf (nnb, rmb, nnb, rvb, rvb) 
call  mply (tm, rvb, rv,nn,  l,nnb) 

c  **  write  rv  to  data  file  (fcoef ficients) 

write (11, 100)  (rv (i) , i=l, nn) 

return 

end 

0  ********************************************************* 

SUBROUTINE  GQUAD (POINT, H) 

c  ** 

c  **  This  subroutine  provides  the  integration  points  and 

c  **  weights  for  the  Gaussian  Quadrature  procedure, 

c  ** 

c  **  variables: 
c  **  ouput: 

c  **  point  —  location  for  numerical  integration 

c  **  weight  —  weight  for  each  point 

c  ** 

0  *********************************************************** 

REAL* 8  H (10) ,  POINT (10) 

C  **  INTEGRATION  POINTS  AND  WEIGHTS  *** 

POINT (1)=. 973906528517172 
POINT (2) =.865063366688985 
POINT (3) =.679409568299024 
POINT (4 >=.433395394129247 
POINT (5) =.148874338981631 
POINT (6) =-.973906528517172 
POINT (7) =-.865063366688985 
POINT (8) =-.679409568299024 
POINT (9)=-. 433395394129247 
POINT (10 >=-.148874338981631 
H(l)=. 066671344308688 
H<2)=. 149451349150581 
H (3) =.219086362515982 
H(4)=. 269266719309996 
H(5)=. 295524224714753 
H(6)=. 066671344308688 
H(7)=. 149451349150581 
H(8)=. 219086362515982 
H(9)=. 269266719309996 
H (10) =.295524224714753 
RETURN 
END 

0  ********************************************************* 
subroutine  bmatrix <a,b, n, zeta,eta,bm) 

c  *  * 

c  **  This  subroutine  forms  the  bmatrix.  The  bmatrix  consists 
c  **  of  three  row  vectors  { [H,yy] / [H,xx] /[H,xy] ) . 
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variables: 
input : 

‘  a  —  plate  length  (in) 

b  —  plate  width  (in) 

'  zeta  —  nondimensionalized  x-location 

:  eta  —  nondimensionalized  y-location 

:  output: 

'  bm  —  bmatrix  of  order  3xN  N= (n+1 ) * (n+1 ) 

r********************************************************** 

implicit  real*8 (a-h, o-z) 
dimension  bm(3, (n+1) * (n+1) ) 

r  initialize  constants  ***************** 

rl=4/b**2 
r2=4/a**2 
r3=4/ (a*b) 
icol=0 


begin  loop  to  formulate  the  bmatrix  *** 

do  10  i=0,n 
do  10  j=0,n 
icol=icol+l 
iyy=j*  ( j-D 

ixx=i* (i-1 ) 
ixy=i*  j 

if  (iyy  .eq.  0)  then 
bm (1, icol) =0 . 0 
else 

bm(l, icol) =iyy*eta** ( j-2) *zeta**i*rl 
endif 

if  (ixx  .eq.  0)  then 
bm<2, icol) =0 . 0 
else 

bm(2, icol) =ixx*zeta** (i-2) *eta**  j*r2 
endif 

if  (ixy  .eq.  0)  then 
bm(3, icol) =0 . 0 
else 

bm(3, icol) =ixy*zeta** (i-1) *eta** ( j —1 ) *r3 
endif 
continue 
return 
end 


subroutine  astif f (a, zeta, tk, theta, nl, break, tcoef , ip, gfac, 
+  ainv,aintv) 

t 

*  This  subroutine  calculates  the  extensional  stiffness 
'  matrix  for  a  specially  orthotropic,  isotropic  plate. 

t 


variables : 
input : 

theta  --  orientation  angle 
nl  —  number  of  layers 
tk  —  lamina  thickness 


'  output: 

*  ainv 

'  aintv 


inverted  stiffness  matrix 
[ai] * [Nthermal] 

*********************************** 


implicit  real*8 (a-h, o-z) 

dimension  ainv (3, 3) , aintv (3) , therm (3) , break (9) , tcoef (4,9) 
dimension  tk (25) , theta (25) , amat (3,3),q(3,3),f(3) 
pi=3. 141592653589793 


convert  zeta  position  to  x  position  and  evaluate  temperature 
curve  fit  at  the  x  position 


x=. 5*a* (zeta+1. ) 
temp=dcsval (x, 8, break, tcoef ) 
dt=temp-75 . 

c  **  call  property  routine  ** 

call  prop(dt,ip,gfac,el,e2, al,a2, rnul2,gl2) 
print*, dt, ip, gfac, el, e2, al, a2, rnul2,gl2 

c  **  initialize  stiffness  matrix  and  thermal  load  vectors 

do  10  i=l, 3 
f (i)=0.0 
do  10  j=l, 3 
amat (i,  j) =0 . 0 
10  continue 

c  **  form  stiffness  matrices 

do  50  ic=l,nl 
c=dcos (theta (ic) ) 
c2=c*c 
c3=c2*c 
c4=c2*c2 

s=dsin (theta (ic) ) 

s2=s*s 

s3=s2*s 

s4=s2*s2 

c  **  compute  reduced  stiffnesses  ** 

rnu21= (e2/el) *rnul2 
qll=el/ (I~rnul2*rnu21) 
ql2=rnul2*e2/ (l-rnul2*rnu21 ) 
q22=e2/ (l-rnul2*rnu21) 
q66=gl2 

c  **  compute  transformed  reduced  stiffnesses  ** 
q(l, 1) =qll*c4+2* (ql2+2*q66) *s2*c2+q22*s4 
q(l, 2) = (qll+q22~4*q66) *s2*c2+ql2* (s4+c4) 
q(2,2)=qll*s4+2* (ql2+2*q66) *s2*c2+q22*c4 
q(l, 3) = (qll-ql2-2*q66) *S*c3+ (ql2-q22+2*q66) *s3*c 
q(2, 3) = (qll-ql2-2*q66) *s3*C+ (ql2-q22+2*q66) *s*c3 
q(3, 3) = (qll+q22-2*ql2-2*q66) *s2*c2+q66* (s4+c4) 

<3(2,1)  =q  (1,2) 
q (3, 1 ) =q (1,3) 
q(3, 2) =q(2, 3) 

c  **  compute  transformed  thermal  expansion  coefficients  ** 
ax=al*c2+a2*s2 
ay=al*s2+a2*c2 
axy=2*s*c* (al-a2) 

c  **  compute  terms  for  thermal  forces  ** 

therm (1) = (q (1, 1) *ax+q(l , 2) *ay+q (1, 3) *axy) *tk (ic) *dt 
therm(2)=*(q(2, 1)  *ax+q(2,2)  *ay+q(2,  3)  *axy)  *tk(ic)  *dt 
therm (3 )  «*  (q(3, 1) *ax+q (3, 2) *ay+q (3, 3) *axy) *tk (ic) *dt 

c  **  form  stiffness  matrix  and  thermal  load  vector  ntv  ** 
do  20  i=l , 3 
f (i) =f (i) +therm(i) 
do  20  j=l, 3 

amat (i, j) =amat (i,j)+q(i, j)*tk(ic) 

20  continue 

50  continue 

do  15  1=1,3 

15  print*,  (amat (i, j) ,  j=l,  3) 

print*, (f (i) , i=l, 3) 

c  **  invert  the  stiffness  matrix  ** 

call  inv (amat, ainv) 

c  **  calculate  ai*f  ** 

call  mply (ainv, f , aintv, 3, 1, 3) 

return 

end 


■k  ★ 
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Q  ★★**★*★★★*★★**★*★★★★*★★*★#★*****★★***★*★★*******★★*★***★★★**•* 

subroutine  prop (at, ip, gfac,el,e2,al,a2, rnul2, gl2) 

q  it  it 

c  **  This  subroutine  accounts  for  the  material  property 
c  **  variation  wrsp  to  temperature, 
c  ** 

c  **  variables: 

c  **  input: 

c  **  dt  —  change  in  temperature 

C  **  ip  —  material  property  specification 

c  **  gfac  —  material  property  variation  with 

c  **  temperature  factor 

c  ** 

c  **  output: 

c  **  el  —  modulus  of  elasticity  (fiber  direction) 

c  **  e2  —  modulus  of  elasticity  (perp.  to  1) 

c  **  al  —  coefficient  of  thermal  expansion  (1) 

c  **  a2  —  coefficient  of  thermal  expansion  (2) 

c  **  rnul2  —  Poisson's  ratio 

c  **  gl2  —  shear  modulus 

Q  ★★★★★★★♦★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★•AT********* 

implicit  real*8 (a-h, o-z) 

c  **  specify  which  material  properties  to  use  ** 
c  **  properties  taken  from  Jones  'Mechanics  of  Composite  ** 
c  **  Materials'  p.70  &  199. 

if  (ip  .eq.  1)  then 

c  **  Graphite/Epoxy  material  properties  ** 

elrt=30 . 0e6 
e2rt“. 75e6 
rnul2=. 25 
gl2rt=. 37be6 
alrt=-0 . 21e-6 
a2rt=16 . Oe-6 
f ac=l . 0-gf ac*dt 
f ac2=l . 0+gf ac*dt 
el=elrt*fac 
e2=e2rt*fac 
al=alrt*fac2 
a2=a2rt*fac2 
gl2=gl2rt*fac 

endif 

if  (ip  .eq.  2)  then 

c  **  Boron/Epoxy  material  properties  ** 

elrt=30 . 0e6 
e2rt=3 . 0e6 
rnul2=. 30 
gl2rt=l . 0e6 
alrt=3 . 5e-6 
a2rt=ll . 4e-6 
f ac*l . 0-gf ac*dt 
f ac2=l . 0+gf ac*dt 
el=elrt*fac 
e2=e2rt*f ac 
al=alrt*fac2 
a2=a2rt*fac2 
gl2=gl2rt*fac 

endif 

if  (ip  .eq.  3)  then 

c  **  E-Glass/Epoxy  material  properties  ** 
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elrt=7 . 8e6 
e2rt=2 . 6e6 
rnul2= . 25 
gl2rt=l . 25e6 
alrt=3 . 5e-6 
a2rt=ll . 4e-6 
fac=l . O-gf ac*dt 
fac2=l . O+gf ac*dt 
el=elrt*fac 
e2=e2rt*f ac 
al=alrt*fac2 
a2=a2rt*fac2 
gl2=gl2rt*f ac 

endif 

if  (ip  .eq.  4)  then 

c  **  Aluminum  material  properties  (isotropic  case)  ** 

elrt=10 . 3e6 
rnul2=. 33 

gl2rt=elrt/ (2* (l+rnul2)  ) 
alrt=12 . 8e-6 
f ac=l . 0-gfac*dt 
f ac2=l . O+gf ac*dt 
el=elrt*f ac 
e2=el 

al=alrt*fac2 

a2=al 

gl2=gl2rt*f ac 

endif 

return 

end 

0  ***★*★★★***★*****★****★★***********★*★★*★★★★★***★★★**★★★★★*** 
SUBROUTINE  INV(XM,XINV) 

c  ** 

c  **  Thi3  subroutine  computes  the  inverse  of  a  3x3  matrix, 
c  ** 

0  A************************************************************* 

IMPLICIT  REAL*8 (A-H,0-Z) 

DIMENSION  XM (3,3) ,XINV(3,3) ,C(3,3) 

C  **  XI=ADJ (X) /DET (X)  ** 

C  **  ADJ (X) “TRANSPOSE  OF  COFACTOR  MATRIX  ** 

C(1,1)=XM(2,2) *XM (3, 3) -XM (2,3)*  XM (3,2) 

c  **  account  for  -  when  taking  the  det  of  x  by  switching  order 

c  **  of  multiplication  for  the  odd  sum  i+j  of  x(i,j) 

C (2, 1 ) =XM (1,3) *XM (3, 2) -XM (1,2) *XM (3,3) 

C ( 3 , 1 ) =XM (1,2) *XM(2, 3) -XM(1, 3) *XM(2,2) 

C (1, 2) =XM(2, 3) *XM(3, 1) -XM (2,1) *XM (3, 3) 

C ( 2 , 2 ) =XM ( 1 , 1) *XM(3,  3) -XM ( 1 , 3) *XM(3,1) 

C (3, 2) =XM (1, 3) *XM (2, 1 ) -XM (1,1) *XM (2, 3) 

C (1 , 3) =XM (2, 1 ) *XM (3,2) -XM (2,2) *XM (3,1) 

C (2, 3) =XM (1,2) *XM ( 3/ 1 ) -XM (1,1) *XM (3,2) 

C(3, 3) =XM(1, 1) *XM(2, 2) -XM (1,2)*  XM (2,1) 

DET=XM (1,1)*C(1,1) +XM (2 , 1 ) *C (2 , 1 ) +XM (3,1) *C(3, 1) 

DO  10  1=1,3 
DO  10  J=l, 3 
XINV(I, J)=C(J, I) /DET 
10  CONTINUE 

RETURN 

END 


<S) 


q  *************************** ********************************** 

subroutine  sdistribution (a, b, n, nx, ny, fcoef , bgrid) 

c  *  * 

c  **  This  subroutine  determines  the  stress  distribution  for 
c  **  the  requested  grid  parameters, 
c  *  * 

c  **  variables: 
c  **  input: 

c  **  a  —  plate  length 

c  **  b  —  plate  width 

c  **  n  —  order  of  polynomial 

c  **  nx  —  #  of  grid  points  in  x-direction 

c  **  ny  —  #  of  grid  points  in  y-direction 

c  **  fcoeff  —  stress  function  coefficients 

c  ** 

c  **  output: 

c  ** 

£  ************************************************************** 
implicit  real*8 (a-h, o-z) 
dimension  bgrid{3, (n+1) * (n+1) ) 
dimension  fcoef ( (n+1) * (n+1) ) ,eload(3) 
nn= (n+1 ) *  (n+1 ) 

c  **  put  origin  of  the  grid  in  the  center  of  the  plate  to  utilize 
c  **  the  symmetry  of  the  problem 

c  dx=a/ (2* (nx-1) ) 

c  dy=b/ (2* (ny-1) ) 

c  xloc=a/2-dx 

c  yloc=b/2 

c  **  to  verify  symmetry  put  the  origin  of  the  grid  at  the  left  edge 

dx=a/ (nx-1) 
dy=b/ (ny-1) 
xloc=0 . 0-dx 
yloc=0 . 0 

c  **  calculate  the  load  intensities  (N)  for  each  grid  point 

do  10  i=l,nx 
xloc=xloc+dx 
zeta= (2*xloc/a) -1 
yloc=0 . 0 
do  10  j=l,ny 
eta= (2*yloc/b) -1 

call  bmatrix (a,b, n, zet a, eta, bgrid) 
call  mply (bgrid, fcoef, eload, 3, 1, nn) 
write  (12, 100)  xloc,yloc,  (eload (ii) , ii=l, 3) 
yloc=yloc+dy 
10  continue 

100  format (lx, f 7 . 3, f 7 . 3, 2x, 3f 18 . 6) 

return 

end 

Q  ***★★*★★***★*****★*★******★*•»  ******************************** 

SUBROUTINE  BTAB (A, B, R, M, N) 

c  ** 

c  **  Thi3  subroutine  calculates  [B] trans* [a] * [ B )  matrix, 
c  *  * 

c  **  variables: 
c  **  input: 

c  **  A  --  matrix  of  order  M*M 

c  **  B  —  matrix  of  order  M*N 

c  **  N  —  order  (n+1)*  (n+1) 

c  **  output: 


c 

c 

c 


*  * 
*  * 


R  —  matrix  of  order  N*N 


*************************************************************** 
implicit  real*8 (a-h, 0-2) 

DIMENSION  A <M,M) ,  B  (M,  N)  ,  R  (N,  N) 

DO  40  1=1, N 
DO  30  J=1,N 
DY=0 . 0 
DO  20  K=1 , M 

IF  (B (K, I)  .EQ.  0.0)  GO  TO  20 
CY=0 . 0 
DO  10  L=1,M 

IF  (B (L, J)  .EQ.  0.0)  GO  TO  10 
CY=CY+A(K,L) *B(L, J) 

10  CONTINUE 

DY=DY+CY*B (K, I) 

20  CONTINUE 

R(I,  J)=DY 
30  CONTINUE 

40  CONTINUE 
RETURN 
END 

c  ************************************************************** 
SUBROUTINE  BTA (L, M, N, A, B, C) 

c  ** 

c  **  This  subroutine  calculates  [B] transpose* [A] . 

c  ** 

c  **  variables: 

c  **  input: 

c  **  A  —  matrix  of  order  N*L 

c  **  B  —  matrix  of  order  N*M 

c  **  output: 

c  **  C  —  matrix  of  M*L 

c  ************************************************************** 
INTEGER  I,J,K,L,M,N 
REAL* 8  A (N, L) , B (N, M) , C (M, L) ,DY 

DO  10  1=1 , M 
DO  10  J=1 , L 
DY=0 . 0 
DO  20  K=1,N 

20  DY=DY+B (K, I ) *A (K, J) 

C(I,  J)=DY 
10  CONTINUE 
RETURN 
END 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  MPLY (A, B, R, M, N, L) 

*  * 

**  This  subroutine  calculates  [aj*[b]  and  stores  it  in  [r] . 
** 


** 
** 
** 
*  * 
*  * 
★  ★ 


variables : 
input : 

A  —  matrix  of  order  M*L 
B  —  matrix  of  order  L*N 
output : 

R  —  matrix  of  order  M*N 


rv - uiattiA  u  i.  i.  n  w 


IMPLICIT  REAL* 8 (A-H, O-Y) 
DIMENSION  A (M, L) ,  B (L, N) ,  R (M, N) 


DO  20  1=1, M 
DO  20  J=1,N 
Y=0.0 

DO  10  K=1 , L 

10  Y=Y+A (I,K) *B(K,  J) 

R  (I,  J)  =Y 
20  CONTINUE 
RETURN 


dll 


END 


0  ************************************************************ 
SUBROUTINE  TRANSFORM (N, TM) 

c  ** 

c  **  This  subroutine  calls  the  transformation  matrix  subroutine, 
c  ** 

c  **  variables: 

c  **  input: 

c  **  N  —  order  of  stress  function  polynomial 

c  **  output: 

c  **  TM  —  transformation  matrix 

q  ************************************************************** 

IMPLICIT  REAL*8 (A-H,0-Z) 

DIMENSION  TM< (N+l) * (N+l) , (N-3) * (N-3) ) 

11=0 

DO  10  1=0, N 
DO  10  J=0, N 
11=11+1 


JJ=0 

DO  10  K=4 , N 
DO  10  L=4 , N 
JJ-JJ+1 

CALL  TMAT (I,K,AFAC) 

CALL  TMAT ( J, L, BFAC) 

TM ( 1 1 , J J ) =AFAC *  BFAC 

10  CONTINUE 
RETURN 
END 

0  ************************************************************* 
SUBROUTINE  TMAT (II, 12, FAC) 

c  ** 

c  **  This  subroutine  forms  the  linear  transformation  matrix 
c  **  to  account  for  the  zero  edge  boundary  conditions, 
c  ** 

c  **  variables: 
c  **  input: 

c  **  II  — 

c  **  12  — 

c  **  output: 

c  **  FAC  — 

0  ************************************************************ 
IMPLICIT  REAL*8 (A-H, O-Z) 

ICOUNT=Il+l 

GOTO (10,20,30,40) ICOUNT 

IF  (II .EQ. 12)  THEN 

FAC=1 . 0 

ELSE 

FAC=0 . 0 

END  IF 

RETURN 

10  N2=(-1)**I2 

N4=<-1) ** (12-1) 

FAC=<-2* (1+N2J+I2* (1-N4) ) /4. 

RETURN 

20  N2=(-1)**I2 

N4= (-1) ** (12-1) 

FAC= (3* (N2-l)+I2* (1+N4 ) ) /4. 

RETURN 

30  N4= (-1) ** (12-1) 

FAC-12* (N4-1) /4 
RETURN 

40  N2»(-1)**I2 

N4=(-l) ** (12-1) 

FAC= ( 1 -N2 -12* (1+N4) ) /4 

RETURN 

END 


d!2 


